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Abstract 

High-density lipoprotein (HDL) is believed to play an important role in lowering cardiovascular disease (CVD) risk by 
mediating the process of reverse cholesterol transport (RCT). Via RCT, excess cholesterol from peripheral tissues is carried 
back to the liver and hence should lead to the reduction of atherosclerotic plaques. The recent failures of HDL-cholesterol 
(HDL-C) raising therapies have initiated a re-examination of the link between CVD risk and the rate of RCT, and have brought 
into question whether all target modulations that raise HDL-C would be atheroprotective. To help address these issues, a 
novel in-silico model has been built to incorporate modern concepts of HDL biology, including: the geometric structure of 
HDL linking the core radius with the number of ApoA-l molecules on it, and the regeneration of lipid-poor ApoA-l from 
spherical HDL due to remodeling processes. The ODE model has been calibrated using data from the literature and 
validated by simulating additional experiments not used in the calibration. Using a virtual population, we show that the 
model provides possible explanations for a number of well-known relationships in cholesterol metabolism, including the 
epidemiological relationship between HDL-C and CVD risk and the correlations between some HDL-related lipoprotein 
markers. In particular, the model has been used to explore two HDL-C raising target modulations, Cholesteryl Ester Transfer 
Protein (CETP) inhibition and ATP-binding cassette transporter member 1 (ABCA1) up-regulation. It predicts that while CETP 
inhibition would not result in an increased RCT rate, ABCA1 up-regulation should increase both HDL-C and RCT rate. 
Furthermore, the model predicts the two target modulations result in distinct changes in the lipoprotein measures. Finally, 
the model also allows for an evaluation of two candidate biomarkers for in-vivo whole-body ABCA1 activity: the absolute 
concentration and the % lipid-poor ApoA-l. These findings illustrate the potential utility of the model in drug development. 

Citation: Lu J, Hubner K, Nanjee MN, Brinton EA, Mazer NA (2014) An In-Silico Model of Lipoprotein Metabolism and Kinetics for the Evaluation of Targets and 
Biomarkers in the Reverse Cholesterol Transport Pathway. PLoS Comput Biol 10(3): e1003509. doi:10.1371/journal.pcbi. 1003509 

Editor: Daniel A. Beard, University of Michigan, United States of America 

Received September 2, 2013; Accepted January 22, 2014; Published March 13, 2014 

Copyright: © 2014 Lu et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted 
use, distribution, and reproduction in any medium, provided the original author and source are credited. 

Funding: This work was supported by the Roche Postdoctoral Fellowship fund of F. Hoffmann-La Roche AG, Basel, Switzerland, project number RPF-153. The 
funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. 

Competing Interests: I have read the journal's policy and have the following conflicts: JL and NAM are employees of F. Hoffmann-La Roche; JL and NAM have 
stock ownership of F. Hoffmann-La Roche. KH, MNN and EAB declare no competing interest. 

* E-mail: james.lu.jl1@roche.com 



Introduction 

Epidemiological studies have shown that high levels of low- 
density lipoprotein cholesterol (LDL-C) as well as low levels of 
high-density lipoprotein cholesterol (HDL-C) are associated with 
increased cardiovascular disease (CVD) risk [1,2]. While LDL-C 
lowering therapies have been shown consistendy to reduce CVD 
risk, there is significant residual risk that remains to be managed 
[2] . The strong inverse association between HDL-C and CVD risk 
has led to the "HDL-C hypothesis", whereby all HDL-C raising 
therapies should be anti-atherogenic [2,3]. Currendy, the anti- 
atherogenic activity of HDL is mainly attributed to its role in 
mediating reverse cholesterol transport (RCT), whereby choles- 
terol is effluxed from peripheral tissues and transported to the liver 
for biliary excretion [4] . However, the recent failures of a number 
of HDL-C raising intervention trials [5-7] have called for a re- 
examination of the HDL-C hypothesis. It has long been thought 
that HDL-C is a reliable biomarker for cholesterol efflux from 
tissues [8]. However, the several recent failed HDL-C raising 



intervention trials provide mounting evidence that at least under 
certain conditions, the plasma concentration of HDL-C, a very 
simple and static measure, is inadequate for characterizing the rate 
of RCT, which is a complex and dynamic process [8] . A revision 
of the HDL-C hypothesis to the "HDL flux hypothesis" has been 
proposed, whereby interventions should be aimed at promoting 
cholesterol efflux to HDL, and hence the overall RCT rate, 
independendy of their effects on HDL-C levels [9,10]. Hence, 
there is now a pressing need to better understand the role of HDL- 
C raising targets in the context of RCT and to identify biomarkers 
which could provide information on the flux rate through the 
RCT pathway [8]. Our modeling effort is focused on addressing 
these issues. 

A number of previous mathematical models have focused on 
various aspects of lipid metabolism; see [11,12] for recent reviews. 
Of the existing models, some describe metabolic processes at a 
mechanistic level [13-19], while others have been empirically 
derived from tracer kinetic studies [20-22]. In general these 
models were built to describe the dynamics of HDL and the other 
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Author Summary 

Epidemiological studies have shown a strong inverse 
association between HDL-C and cardiovascular risk and 
led to the formulation of the "HDL cholesterol hypothesis": 
under this hypothesis, interventions raising HDL-C should 
decrease risk. However, the recent failures of HDL-C raising 
therapies in improving cardiovascular disease risk in 
outcomes trials have suggested a need to revise the 
hypothesis to account for the contrary data. An "HDL flux 
hypothesis" has emerged: it is not HDL-C level per se which 
forms the basis for reducing risk, but it is the flux rate of 
reverse cholesterol transport that drives risk reduction. We 
propose that, the concentration of HDL cholesteryl ester in 
plasma simply reflects the ratio of input rate of reverse 
cholesterol transport into the HDL compartments to its 
clearance rate. A challenge in identifying targets under the 
new conceptual framework is the feedback process that 
occurs between the input rate and the clearance rate of 
HDL-C. To meet this challenge, we have built a systems 
model which incorporates the main processes of HDL 
metabolism to elucidate the relationships between target 
modulations and the reverse cholesterol transport rate. 



major lipoprotein classes, which include LDL, intermediate 
density lipoprotein (IDL) and very low density lipoprotein (VLDL), 
describing lipid transport between these particles mediated by the 
cholesteryl ester transport protein (CETP) in the normal or basal 
state, and the effects of genetic mutations and/ or drug interven- 
tions on these processes. While valuable insights have been gained 
from these models, none can be used to predict the associated 
changes in the RCT rate since they lack a mechanistic description 
of ApoA-I dynamics and other key processes involved in the RCT 
pathway. The latter include the lipidation of lipid-poor ApoA-I via 
its interaction with ATP-binding cassette transporter member 1 
(ABCA1), the key process in the initiation of RCT [8], as well as 
processes of HDL remodeling which lead to the delivery of 
cholesterol from HDL to other lipoproteins and cells, and the 
regeneration of lipid poor ApoA-I [23]. In all the existing models 
except the ones by Htlbner et al [16] and Adiels et al [24], the 
dynamics of apolipoproteins that cover the surface of lipoprotein 
particles are not described. While each VLDL, IDL and LDL 
particle contains only one ApoB molecule per particle, for HDL 
particles the number of ApoA-I molecules per particle may vary 
from 2 to 4 or more depending on HDL size [25]. This variation 
results from HDL remodeling processes such as particle fusion, 
CETP-mediated lipid transport, lipolysis and esterification where- 
by particles can gain or lose core lipid content as well as ApoA-I 
molecules [23]. While it has been shown experimentally and 
theoretically that the number of ApoA-I molecules on a given 
HDL particle is intrinsically linked to the particle size [25], this 
important relationship has yet to be incorporated into a 
mechanistic model of HDL metabolism. 

In this paper, we propose a novel model of lipoprotein 
metabolism and kinetics (the LMK model) that provides an 
integrated description of the dynamics of cholesterol and ApoA-I 
in plasma. In particular, the model captures the initiation of RCT 
from the lipidation of lipid-poor ApoA-I by the ABCA1 
transporter, the generation of nascent discoidal and nascent 
spherical particles, HDL particle fusion, CETP mediated lipid 
transfer between HDL and other lipoproteins, and the dissociation 
of excess ApoA-I from mature spherical a-HDL due to remodeling 
processes. The model is calibrated to: lipoprotein measures for 
normal and CETP deficient subjects; cholesteryl ester (CE) and 



ApoA-I fluxes measured in normal subjects; data on the fractional 
catabolic rate (FCR) of ApoA-I. The structure and the kinetic 
constants of our model provide an explanation for the relationship 
between FCR of ApoA-I and HDL particle size. To our 
knowledge the LMK model is the first to provide a mechanistic 
basis for the linkage between the metabolism of ApoA-I and the 
cholesterol component of HDL. The model has been validated by 
simulating patients with genetic mutations in the HDL metabolism 
pathway and the predictions are compared with lipoprotein 
measures reported in literature. Finally, the model was used to 
evaluate targets that could potentially increase RCT and to 
identify relevant biomarkers, as part of the effort to support drug 
discovery and development using a model-based approach. 

Results/Discussion 

Model structure 

The LMK model is shown schematically in Figure 1 , focused on 
the RCT pathway and a number of targets contained within it, for 
instance CETP, ABCA1, ApoA-I and SRB1. The LMK model 
describes the synthesis of ApoA-I and the initiation of RCT by the 
interaction of lipid-poor ApoA-I with ABCA1 leading to the 
formation of mature, spherical a-HDL. The HDL remodeling 
processes represented in the model include: the fusion of spherical 
HDL particles (arrow 5 of Figure 1); the exchange and elimination 
of CE in spherical HDL by interaction with CETP (arrows 12-14) 
and SRB 1 (arrow 7); the regeneration of lipid-poor ApoA-I from 
spherical HDL particles (arrow 3). Lipid-poor ApoA-I is assumed 
to be eliminated via the kidney (arrow 4), while the spherical HDL 
particles are assumed to be eliminated by a holo-uptake 
mechanism with a rate dependent on the particle size (arrow 6). 
The transfer and elimination of CE in LDL and VLDL pools are 
also represented (arrows 9-11). Our approach is to adequately 
describe the metabolic processes, while keeping the model as 
simple as possible. The representations of lipoprotein components 
and metabolic processes in the LMK model reflect these 
principles. 

Lipoprotein representation. While HDL particles are 
heterogeneous in size and composition [8], for the purposes of 
understanding RCT we only consider two HDL particle classes: 
spherical, a-HDL and small, lipid-poor ApoA-I. Amongst the 
apolipoproteins and lipid species contained in 0£-HDL particles, 
the LMK model has an explicit representation of ApoA-I and CE. 
Although there are a large number of species in the HDL 
proteome (e.g., ApoA-II, ApoE) and lipidome (e.g., triglycerides, 
phospholipids) which may be relevant in particular diseased states, 
they play a secondary role in the characterization of RCT. We 
make the assumption that the protein moiety contains 60% ApoA- 
I by weight, with all other proteins contributing the remaining 
40%. This is within the range of values reported in literature 
[25,26]. Under this assumption, the total concentration of ApoA-I 
is represented as an explicit variable that changes as a direct result 
of the metabolic processes described in the model while ApoA-II 
and other HDL apolipoproteins are implicit quantities: namely, 
they are assumed to change in concert with ApoA-I so as to keep 
the weight fraction constant. Similarly, the LMK model explicidy 
represents CE in the particle classes of a-HDL, VLDL and LDL 
but represents TG in a-HDL only implicitly. That is, the ratio of 
TG/CE in a-HDL particles is assumed to be 13% which is 
consistent with the range of values reported in healthy subjects 
[25,27]. The amounts of free cholesterol (FC) and phospholipids 
(PL) per HDL particle are implicitly represented in the LMK 
model: they depend on the a-HDL size, in a manner analogous to 
the treatment of PL in [16]. In particular, given the CE content of 
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Figure 1 . A schematic representation of the model. The arrows shown in the diagram denote the processes represented by the model and the 
boxes with italicized text denote mediators that are explicitly represented. The process arrows are numbered, refering to the reaction number shown 
in Table 2. The arrows leading from the nascent sphere towards the a-HDL pool represent the 2 scenarios that may occur in the transformation of 
newly formed particles: they may either enter the a-HDL pool as distinct particles (the dashed arrow) or fuse with the existing ones (solid arrow). 
doi:10.1371/journal.pcbi.1003509.g001 



an a-HDL particle its core size can be inferred and the FC and PL 
content on the surface can be computed using the updated Shen 
model; see [25]. With our choice of lipoprotein representation, the 
species represented in the model are given in Table 1. 

Metabolic processes. The full list of reactions represented in 
the LMK model (as schematized in Figure 1) is shown in Table 2. 
We would like to point out that the remodeling flux (arrow 3 of 
Figure 1) based on geometric concepts developed in [25] is an 
original contribution of our work. The remodeling flux expression, 
^ i rem(CE ( j(f),^4i X (i),JV a (f)), represents the excess ApoA-I within the 
pool of a-HDL particles given the core cholesteryl ester content, 
particle concentration and the amount of ApoA-I covering the 



Table 1. Species represented in the model. 





Symbol 


Description 


Units 


A\ P 


Lipid-poor ApoA-I 


mg/dL 


A* 


ApoA-I in the a-HDL pool 


mg/dL 


N a 


Particle concentration of cc-HDL 


mmol/dL 


CE a 


cholesteryl ester in a-HDL 


mg/dL 


CE LDL 


cholesteryl ester in LDL 


mg/dL 


CEvLDL 


cholesteryl ester in VLDL 


mg/dL 



doi:1 0.1 371 /journal.pcbi.1 003509.t001 



surface. Its derivation based on geometric concepts of a-HDL 
particles is discussed in more detail in the Methods section. The 
holo-uptake of a-HDL particles is thought to be mediated by a 
number of receptors, which are not well understood [28,29]. In 
order to account for the possible size-dependence in the uptake 
rate of a-HDL particles, the functional dependence kh 0 i o (d) is 
utilized. This is also discussed in more detail in the Methods 
section. 

The model constants are shown in Table 3 while the list of 
parameters are given in Table 4; the prior values of parameters 
and their posterior estimates are discussed in the next section. 
With the list of reactions given in Table 2, the LMK model can be 
expressed as the following system of ODEs: 

dA, p (t) = 

dt ~ (la) 
r l ?-kABCA\Ai p (t)-k kidas . y Ai p (t) +k dissoc F tem (CE a (t),A ol (t),N lt (t)) 



dA a (t) = 

dt (lb) 

^ABCAl^/ P (0-^dissoc-F'rem(CE a (0,^a(0^a(0)-^holo(rfM ;I (0 
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Table 2. Reactions represented in the model. 





# 


Reaction 




Description 


Rate expression 


Ref. 


1 






ApoA-l synthesis 


fi 

in 


[22,88] 


2 


A, ->A +v mc C¥ 4 


2m A 


Initiation of RCT by interaction with ABCA1 




[79,80,89] 




4 —* A, 




Regeneration of lipid~poor ApoA - l via HDL remodelincj 




[11 1^ 7fi Qfl Q91 


4 






Kidney removal of lipid-poor ApoA-l 


^kidney A/p(0 


[43,48] 


5 






Fusion of nascent spherical particles with mature a-HDL 


^ABCAlAW kfN,(f) 
2m A l+k f N,(t) 


[23,93] 


6 


{A„CE„N,}->0 




HDL particle holo-uptake 


A "holoW x K(0.CE.(0,iV„(f)} 


[28,29] 


7 


CE.-.0 




SR-B1 mediated removal of CE from HDL particles 




[4,8,18,36] 


8 


0^CE VL DL 




Synthesis of CE in VLDL 


( .VLDL 
in 


[20,21] 


9 


CEyLDL ~* 0 




Elimination of CE from VLDL 


/f out Lt VLDL ( '> 


[20,21] 


10 


CE LDL ^0 




Elimination CE from LDL 


fc oS LcE LD L W 


[20,21] 


11 


CEyLDL ^CE LDL 




VLDL conversion to LDL via lipolysis 


/f VL CE VLDL ( " 


[20] 


12 


CE a ^-CE V LDL 




CETP mediated CE transfer from HDL to VLDL 


A-£ET p CE lW 


[13,17,18,20,30] 


13 


CE a ->CELDL 




CETP mediated CE transfer from HDL to LDL 


feg ETP CE«(0 


[13,20,30] 


14 


CEldl-^CE^ 




CETP mediated CE transfer from LDL to HDL 


/f LH TPcE LDL« 


[13,18,20,30] 



doi:1 0.1 371 /journal.pcbi.1 003509.t002 



dNM 
dt 



1 



2m A l+k f N x (t) 



A, p (t)-k ho \ 0 (d)N a {t) (lc) 



dCE<z(t) m c , ... 

— -T. — =y — KABCAl^p(0 + 
at 



/c™ TP CE LDL (?)- [/^v TP +^! TP + ^rb L i +fehoioW]CE a (0 



(Id) 



dt 

*vlCEvldl(0 +/c™ tp CE«(0 ■ 



(le) 



( ,CETP 
v K LH 



+ ^ L ° L )CE LDL (/) 



^CEvldlCQ _ 
dt 

>l LDL + k^TCE^t) - (£ VL +C t DL )CE V LDL(0- 



(If) 



One important quantity that the LMK model can help to assess is 
the rate of reverse cholesterol transport (RCT) at the whole-body level. 
This quantity is thought to play an important role in determining 
cardiovascular disease risk, but is experimentally challenging to 
assess. Using the LMK model, we are able to quantify the flux rate 
of free cholesterol into the nascent disc particles mediated by 
ABCA1: in particular, this is given by 



trie 

RCT rate = plasma volume x y — fcABCAi^fo(0- (2) 

m A 

The term A:abcai^//j(0 represents the transformation rate of lipid- 
poor ApoA-I to nascent discs, which subsequently enter the 



a-HDL pool. The parameter y describes the number of cholesterol 

mc 

molecules per ApoA-I in the nascent discs. converts the 

m A 

molecular mass of ApoA-I to cholesterol. Finally, the volume of 
plasma converts RCT rate to the whole-body level: we assume that 
plasma volume = 3. 15 L in a 70 kg adult [22]. As illustrated in 
Figure 1, nascent discs are transformed into nascent spheres (as 
mediated by the LCAT enzyme [30]) which are assumed to have 
in their cores y CE molecules per ApoA-I. Hence, the RCT 
expression (2) also represents the input rate of HDL-CE into the 
plasma 0£-HDL pool. Note that the factor 2 in the expression for 
reaction 2 (initiation of RCT) accounts for the assumption that 
there are 2 ApoA-I molecules per nascent HDL particle. 

Model calibration 

Parameter estimates: Prior and posterior. The Bayesian 
approach for parameter estimation is a well established method- 
ology which has found applications in various fields of science 
[31], including parameter estimation for models of cellular 
processes [32,33] as well as pharmacokinetics and pharmacody- 
namics (PK/PD) [34,35]. Under this framework, it is assumed that 
a prior distribution is available for (some) parameters as a result of 
previous experimental studies. In combination with calibration 
data, the posterior distribution for the parameters is obtained. 

For most of the LMK model parameters, prior estimates are 
available from literature studies; a detailed discussion of the 
references from which parameter estimates and their uncertainties 
are obtained is given in the Methods section. Using the model 
calibration procedure as discussed in the Methods section, the 
prior is combined with calibration data to give rise to the posterior 
estimates. A list of the prior and posterior values of parameters is 
given in Table 5. It is worthwhile noting that, for the most part, 
the maximum a posteriori (MAP) estimate obtained by the calibration 
process does not depart significandy from the prior. This indicates 
that the calibration data are fairly consistent with the prior 
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estimates. One exception is the parameter ^gRgij which is 
increased significantly from its prior beyond the 1 SD value. This 
result is in agreement with experimental evidence that SRB 1 plays 
a significant role in mediating HDL-CE removal from HDL 
particles [36], in contrast to the expectation of a previous tracer 
kinetics study [21]. The discrepancy may be attributed to the 
limitation of tracer kinetics studies (for instance, [21]) to be able to 
fully identify the SRB1 contribution. Finally, it can be seen that 
the calibration data are sufficiently informative to allow relatively 
precise estimates for the parameters k^ olo and kf, for which there 
was no prior information. It is worth noting the negative sign in 
the estimate for ^olo> which implies that the a-HDL holo-particle 
uptake rate decreases with particle size. The sign of this size- 
dependence is consistent with the hepatic endocytic receptor 
(mitochondrial ATP synthase subunit f!) having a higher affinity 
for the (smaller) HDL-3 as compared to the (larger) HDL-2 
[29,37,38]. 

Calibration data and model explanatory power. In order 
to identify parameter values using the Bayesian approach, 
calibration data are needed. However, the choice of calibration 
data should be made not only for the purpose of quantifying 
parameters, but also with the consideration for the potential utility 
and the explanatory power of the model. More specifically, 
choosing the right types of calibration data can help to increase 
confidence in a model's predictions of specific scenarios; in 
addition, model calibration is also an opportunity to test if the 
model structure, together with prior information on the parameter 
values, can explain important features of the system being studied. 

In our current work, the LMK model was used to explain the 
effects of CETP inhibition on ApoA-I level as well as the inverse 
relationship between the FCR of ApoA-I and particle size. Based 
on this, the calibration data were chosen. In Figures 2, 3 and 4 we 
show the calibration data superimposed with the model simula- 
tion, using the maximum a posteriori parameter set identified by the 
calibration procedure (as described in the Methods section). In 
Figure 2, the decrease in CETP level from 100% to 0% of the 
nominal subject was simulated by decreasing the three parameters 

associated with CETP activity (^HV TP ^HL TP '^LH TP) by the 
same factor. In particular, panels A and B show that the rise in 
HDL-C (the concentration of HDL-C is computed by summing 
CE a and free cholesterol pool, as discussed in the Methods section) 
and ApoA-I in heterozygotes and homozygotes with CETP 
deficiency are fairly well captured by the LMK model; the main 
discrepancy is the under-prediction of HDL-C for CETP 
heterozygotes. To our knowledge, the increase in ApoA-I under 
CETP deficiency or inhibition has not yet been explained by 



Table 3. Model constants. 



Constant 


Description 


Unit 


Value 


Ref. 


m A 


Molecular weight of ApoA-I 


g/mol 


28500 


[16] 


me 


Molecular weight of cholesterol 
(free and esterified) 


g/mol 


386 


[16] 


V CE 


Molecular volume of cholesteryl 
ester 


A 3 


1068 


[25] 


V TG 


Molecular volume of triglyceride 


A 3 


1556 


[25] 


t 


Thickness of HDL surface 


A 


20.2 


[25] 



By convention, the cholesteryl ester mass is measured by quantifying the 
equivalent mass of free cholesterol. 
doi:1 0.1 371/journal.pcbi.1 003509.t003 



existing models: by incorporating the geometric ideas proposed in 
[25] the LMK model provides, for the first time, a way to connect 
the metabolism of ApoA-I and HDL-C. We remark that the LMK 
model is focused on HDL rather than the metabolism of ApoB- 
containing particles, which include LDL and VLDL. In particular, 
the LMK model predicts negligible concentrations of LDL-CE 
and VLDL-CE in CETP homozygotes, which are inconsistent 
with the reported concentrations in these subjects [39-41]. We 
believe that this discrepancy between the LMK model prediction 
and reality is due to /J-LCAT activity [42], which in CETP 
homozygotes could compensate for the lack of CE influx from 
HDL particles by converting free cholesterol on the surface of 
ApoB-containing particles into cholesteryl ester. Finally, Figure 3 
shows that the CE fluxes in the LMK model are consistent with 
the values measured [20], in particular the CE flux from HDL to 
LDL is close to that from LDL to HDL. 

Another important and robust finding that has been observed in 
HDL metabolism is the relationship between FCR of ApoA-I and 
particle size (estimated using a surrogate measure) seen in normal 
subjects [43,44]; in addition, heterozygotes and homozygotes of 
CETP deficiency are also observed to have a decreased FCR of 
ApoA-I [45]. Thus, an important objective of the calibration 
process is to test whether the structure of the model, together with 
an assumption on the linear size-dependence of HDL holo-particle 
uptake rate, can explain this relationship. The inverse relationship 
observed between the FCR of ApoA-I and the ratio HDL-C/ 
ApoA-I (a surrogate measure of HDL size) is shown in Figure 4: in 
particular, data from Schaefer et al [44], Brinton et al [43] and 
Ikewaki et al [45] are given. A linear fit was carried out using the 
pooled data of Schaefer et al [44] , Brinton et al [43] and normal 
subjects from Ikewaki et al [45], with the mean shown as a dashed 
pink line and the 1 SD confidence region shown as dashed black 
lines in Figure 4. The ApoA-I FCR for CETP homozygotes (who 
have large HDL particles) are assumed to be the lowest level 
attainable, hence this value was taken as the "floor" of the fit. The 
LMK model was calibrated to the piecewise linear relationship 
represented by the pink line and Figure 4 shows that simulations 
for normal and CETP mutation subjects (denoted by the asterisk 
symbols) are all in good agreement with the inverse relationship. 
The LMK model reproduces the dependence of ApoA-I FCR on 
CETP primarily by changing the distribution of ApoA-I between 
the lipid-poor and a-HDL pools (which have different clearance 
rates), with a minor contribution from the explicit size dependence 
of holo-particle uptake, k^ 0 \ 0 {d). 

Model validation 

In order to increase confidence in its predictions, the LMK 
model has been validated by simulating a number of scenarios that 
have not been used in the calibration process. In particular, since 
ABCA1 and ApoA-I are important targets in the pathway, the 
literature data on subjects with mutations in these genes [46,47] 
are compared against the model simulations. The heterozygotes 
and homozygotes of ABCA1 mutation are simulated by setting 
^ABCAl (representing ABCA1 activity) to 50% and 0% of its 
nominal value respectively; similarly, heterozygotes and homozy- 
gotes of ApoA-I mutation are simulated by setting the parameter 

rP (representing ApoA-I synthesis rate) to 50% and 0% of its 

nominal value respectively. Figure 5 shows the mean and 95% 
confidence intervals of the model simulations, compared to the 
literature data (mean and SD are given). An examination of the 
results for the heterozygotes shows that, encouragingly, the LMK 
model is able to differentiate between the effects of ABCA1 and 
ApoA-I mutations on HDL-C and ApoA-I levels: both quantities 
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Table 4. Model parameters. 





Description 


Unit 


J? 

m 


Synthesis rate of ApoA-l 


mg/dL/day 




Rate of kidney elimination 


pool/day 


aissoc 


Dissociation rate of excess ApoA-l 


pool/day 


A ABCA1 


Rate constant in the lipidation of lipid-poor ApoA-l via ABCA1 


pool/day 


7 


Stoichiometry of FC to ApoA-l in nascent discs 


unitless 


/.CETP 

K HV 


Rate constant of CE transfer: HDL to VLDL 


pool/day 


,.CETP 

"hl 


Rate constant of CE transfer: HDL to LDL 


pool/day 


,.CETP 
LH 


Rate constant of CE transfer: LDL to HDL 


pool/day 


k VL 


Rate constant of CE transfer: VLDL to LDL 


pool/day 


r VLDL 
in 


Synthesis rate of CE to VLDL 


mg/dL/day 


/.VLDL 
K out 


Rate constant of CE elimination from VLDL 


pool/day 


iLDL 
"out 


Rate constant of CE elimination from LDL 


pool/day 


/.HDL 
A SRB1 


Rate constant of SRBI-mediated CE elimination from HDL 


pool/day 


k c 
holo 


Constant contribution to the rate of a-HDL holo-particle uptake 


pool/day 


^holo 


Size-dependent contribution to the rate of k-HDL holo-particle uptake 


pool/day/nm 


k f 


Parameter governing the particle concentration dependence of fusion rate 


1/(mmol/dL) 



doi:1 0.1 371 /journal.pcbi.1 003509.t004 



decrease more for ApoA-I heterozygotes as compared to ABCA1 
heterozygotes. Furthermore, the LMK model predicts that 
heterozygotes of ABCA1 mutation have smaller HDL particles 
(data not shown), consistent with the data of Asztalos et al [46] . 

Most of the calibration data are static in nature, hence it is of 
particular interest to perform dynamic simulations of the LMK 
model and compare them to existing data. As a validation, we 
would like to see if the LMK model reproduces the characteristic 
biphasic decay curves seen in tracer kinetic experiments with 
labelled ApoA-I. In the LMK model, the injection of radio- 
labelled dose is represented by a small addition to the pool of lipid- 
poor ApoA-I and the fractional dose remaining in the sum of the 
two pools of ApoA-I is plotted; refer to the Methods section for the 
details of the simulation methodology. This is simulated using the 
parameters identified for the nominal subject and the result is 
shown in Figure 6: it can be seen that the simulated decay curve is 
biphasic and similar to the data obtained by digitizing Figure 3 of 
Ikewaki et al [48]. Furthermore, the mean residence time (which is 
the inverse of FCR) of labelled ApoA-I computed from the model 
simulation is 4.2 days, which is in good agreement with the result 
of 4.8±0.3 days as measured in 4 subjects by Ikewaki et al [48]. 

Explaining epidemiological relationship using a virtual 
population 

Having calibrated and validated the LMK model, we use it as a 
platform for exploring the observed epidemiological relationship 
between HDL-C and CVD risk. For this purpose, a virtual 
population is generated in a manner analogous to that of reference 
[49]. In particular, model parameters are sampled from a 
multivariate normal distribution and for each set of parameters 
the "phenotype" of the corresponding virtual subject is simulated 
using the LMK model. As there is no information available on the 
correlation between model parameters in a real population, we 



have assumed them to be uncorrelated and each is drawn from a 
normal distribution with a relative SD = 1 5 % around the value 
corresponding to the posterior values for the nominal subject (see 
Table 5). 

Despite the fact that the parameter distribution in the virtual 
population is uncorrelated, some of the simulation outputs show 
significant correlations as a result of the model structure. Of 
particular interest is the correlation between RCT rate (as defined in 
(2)) and plasma biomarkers. Shown in Figure 7 is the relationship 
between RCT rate and HDL-C within the virtual population: it can 
be seen that there is a surprisingly strong correlation between the two 
quantities (r=0.95). We note that the RCT rate given in (2) 
corresponds to the input rate of HDL-CE into plasma: in fact, the 
plasma concentration of HDL-CE can be expressed as the following: 



HDL — CE = 



RCT rate 



Clearance rate of HDL — CE 



(3) 



where the clearance is defined as the plasma volume multiplied by 
the sum of elimination rate constants. In the LMK model, 
elimination processes for HDL-CE include those mediated by 
CETP and SRB1, as well as the holo-particle uptake. While the 
RCT rate shows a strong correlation with HDL-C, we see that in 
Figure 8 the clearance of HDL-CE shows a much weaker negative 
correlation with HDL-C (r= —0.32). Hence, the simulation results 
suggest that the variation in HDL-C within the virtual population is 
largely attributed to variations in the RCT rate and not due to its 
clearance. Under the "HDL flux hypothesis" [2] that low RCT rate 
results in high CVD risk, the relationship shown in Figure 7 provides 
a plausible explanation for the epidemiological association 
between HDL-C and CVD risk. The same set of virtual subjects is 
also used in subsequent sections for target evaluation and biomarker 
identification. 
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Table 5. Prior and posterior estimates of model parameters 
corresponding to the "nominal subject." 



Parameter 


Prior (mean±SD) 


Posterior (mean±SD) 


in 
ID 

r- 
m 


A 1 .44 ± 1 . 1 o 


ZiJ.4D± 1 . 1 3 


^kidney 


5.19±2.60 


2.42±0.78 


^dissoc 


174±312 


170±191 


ABCAl 


C\C "1 A -+- 1 ~l CC 

yO.Z4_ 1 /.DD 


yD. 1 i3± 1 D./d 




7.55 ±3.94 


10.17±2.19 


A .CETP 
HV 


1.47 ±0.58 


1.49±0.24 


^.CETP 
HL 


5.47±2.05 


6.92±0.81 


.CETP 
LH 


1.98 ±0.70 


2.89±0.34 


k VL 


7.52 ±0.94 


7.70±0.84 


r VLDL 
in 


0.96 ±0.46 


1.50±0.45 


,.VLDL 
out 


0.88±0.37 


1.30±0.35 


tLDL 
K out 


0.67 ±0.08 


0.64±0.07 


/HDL 
"SRB1 


0.31 ±0.12 


0.60±0.08 


^holo 


0.14±0.026 


0.13±0.022 


k 1 
holo 


0 


-0.016±0.004 


kf 


0 


5000±1544 





The prior parameter values were estimated from literature as discussed under 
the Methods section. The SDs of the posterior distribution were estimated using 
the Fisher Information Matrix. 
doi:1 0.1 371/journal.pcbi.1 003509.t005 



HDL-C raising therapies 

The LMK model can be used to evaluate actual and potential 
HDL-C raising therapies, by modulating targets of interest. We 
have used simulated the model for both the nominal subject, as 
well as for a virtual population. 

CETP inhibition. The model predictions of CETP inhibition 
on the nominal subject, together with 95% confidence intervals, 
are shown in Figure 9. In alignment with the calibration data, as 
CETP level decreases both HDL-C and ApoA-I increase strikingly 
(see panels A and B). These changes are associated with a 
significant increase in HDL size as well as a small increase in HDL 
particle concentration (HDL-P) (see panels C and D). Both the 
absolute concentration of lipid-poor ApoA-I and RCT rate remain 
essentially unchanged (see panels E and F). 

To further illustrate the effect of CETP inhibition, in Figure 10 
we simulate the therapy in a virtual population. In particular, we 
select subjects with low HDL-C (40 mg/dL or less) for treatment 
with a hypothetical drug that inhibits the plasma CETP by 80% 
and simulate the changes in HDL-C and RCT rate in the treated 
subjects. It can be seen that the rise in HDL-C does not 
correspond to an increase in RCT rate. In fact, the effects induced 
by CETP inhibition depart from the baseline relationship. This is 
an illustration of a target impacting a biomarker which is 
correlated but not causally linked with the disease mechanism: 
the hypothetical drug does not bring about a therapeutic effect of 
increasing RCT rate despite increasing HDL-C. However, there 
could be a potential CV benefit due to a small (14%) decrease in 
LDL-C (data not shown). 

ABCAl up-regulation. The model predictions for ABCAl 
up-regulation on the nominal subject, together with 95% 



confidence intervals, are shown in Figure 1 1 . The simulation 
results show that as ABCAl activity increases, both HDL-C and 
ApoA-I increases (see panels A and B). Panels C and D of Figure 1 1 
show that these increases reflect not only an increase in HDL size, 
but also increases in particle concentration. In stark contrast to 
CETP inhibition, under ABCAl upregulation the RCT rate is 
predicted to increase markedly (panel E) and the absolute 
concentration of lipid-poor ApoA-I decreases. 

We next consider ABCAl up-regulation for the virtual 
population as shown in Figure 7. In particular, we select the 
same subjects with low HDL-C as was previously chosen for 
CETP inhibition. We simulate a hypothetical drug that increases 
ABCAl activity in each of the treated subjects by 100% and 
examine the changes in HDL-C and RCT rate. As shown in 
Figure 12, under ABCAl up-regulation both HDL-C and RCT 
rate increase. In particular, the changes induced by ABCAl up- 
regulation are predicted to follow the baseline epidemiological 
relationship. 

In order to further elaborate on the differences between the two 
target modulations, we compare in Figure 13 the changes in 
biomarkers for CETP inhibition and ABCAl up-regulation. The 
simulation results show that for a given fold change in HDL-C, 
CETP inhibition gives rise to larger particle sizes but fewer particle 
numbers as compared to ABCAl up-regulation. Due to the 
differences in the particle size and number under the two target 
modulations, for a given fold-change in HDL-C, ApoA-I is 
predicted to increase more under ABCAl up-regulation as 
compared to CETP inhibition. Reassuringly, the simulated 
increases in ApoA-I for CETP inhibition are in fair agreement 
with literature data for the three CETP inhibitors, Dalcetrapib 
[50], Torcetrapib [51] and Anacetrapib [52]. For CETP 
inhibition the predicted decline in LDL-C as the fold-change in 
HDL-C increases is in good agreement with data on the three 
CETP inhibitors (Figure 13, panel E). Conversely for ABCAl up- 
regulation the LMK model predicts an increase in LDL-C with 
increasing fold-change in HDL-C. This finding is a consequence of 
the first-order CETP-mediated transfer processes between HDL, 
VLDL and LDL particles. It is qualitatively consistent with the 
GWAS study of Voight et al [53] which showed that ABCAl SNP 
rs3890182 raised HDL-C and LDL-C by comparable amounts. 

Lipoprotein biomarkers 

A number of studies have shown that CVD risk is correlated 
with plasma biomarkers such as HDL-C [1], HDL-P [54] and pre- 
P\ [55,56] levels. In addition, the combination of NMR analysis of 
HDL with genotyping has also given a glimpse into the possible 
genes associated with HDL particle measures [57]. However, the 
mechanistic basis for these experimental observations as well as 
what underlies the correlations between the plasma biomarkers are 
not well understood. Using the proposed LMK model, we can 
reproduce and explain the correlations between these plasma 
biomarkers. 

In addressing these questions, the simulated biomarkers within 
the population of 2000 virtual patients (as previously shown in 
Figure 7) were studied. The correlation between HDL-P and 
HDL-C within this set of virtual patients (r = 0.74) is shown in 
Figure 1 4, panel A; we see that the simulation result is qualitatively 
similar to the positive correlation shown by Mackey et al [54] (the 
absolute values of HDL-P obtained by NMR are approximately 2- 
fold greater than our simulations which are based on the updated 
Shen model; the discrepancy is discussed in [25]). A positive 
correlation also exists in the virtual population between HDL size 
and HDL-C (r = 0.71), consistent with Mackey et al [54] (see 
Figure 14, panel B). 
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Figure 2. The fit of the model to the calibration data for CETP deficiency: HDL-C (panel A), ApoA-l (panel B), LDL-CE (panel C) and 
VLDL-CE (panel D). The data are as shown in Table 6, obtained by pooling HDL-C and ApoA-l data from references [81-83] and LDL-CE, VLDL-CE 
data from references [20,21]. The model simulation curves were obtained by decreasing the 3 parameters representing CETP activity 
(^hv TP '^hl TP '^lh TP ) f rom 100% to 0% of those corresponding to the nominal subject. 
doi:10.1371/journal.pcbi.1003509.g002 



Due to the growing appreciation for the importance of RCT 
[2,8,58], there are on-going efforts in trying to quantitatively 
assess the steps involved in the process. The ABCA1 transporter 
is involved in the first step of RCT by removing cholesterol from 
peripheral tissues to plasma and its activity level in patients has 
been studied [59]. In particular, ABCA1 gene expression and 
protein concentration on leukocytes has been measured in 
patients with type 2 diabetes, where the data suggested a negative 
correlation between ABCA1 expression and HbAlc levels [59]. 
While there are assays that can quantify ABC A 1 protein levels in 
specific cell types [60], an experimental technique for the 
assessment of ABCA1 activity in-vivo at the whole body level 
has yet to be developed. Given the current experimental 
limitations, there is an interest to evaluate the potential 
effectiveness of plasma-based biomarkers for quantitatively 
assessing ABCA1 activity. 



Using the LMK model, we evaluated the potential effectiveness 
of two biomarkers for ABCA1 activity: firstly, the absolute 
concentration of lipid-poor ApoA-I; secondly, the relative 
concentration of lipid-poor ApoA-I as the percentage of total 
ApoA-I. Figure 15 panel A shows that the former is only weakly 
correlated with ABCA1 activity. In contrast, panel B shows that 
the latter exhibits a strong inverse correlation with ABCA1 
activity; in fact, given a measured value of % lipid-poor ApoA-I, 
the relationship can be used to estimate ABCA1 activity. This 
result can be better understood by the following analysis. From 
equation (la), the absolute concentration of lipid-poor ApoA-I at 
steady state can be expressed as: 

_ ^ + fcdissoc f rem(CE„A,Ay 
/C kidney + /C ABCA1 
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Figure 3. The fit of the model to the calibration data: CE fluxes. The data are as shown in Table 8, taken from reference [20]. The model 

simulation is produced using the point estimate of parameters for the nominal subject. 

doi:10.1371/journal.pcbi.1003509.g003 



C ABCA1 



(4b) 



On the other hand, from (lb) the % lipid-poor ApoA-I can be 
expressed as the following: 



A lp _ fchoioW+fcdjssoc^remCCE^AjAy/A 
A « ^ABCAl 



fe ( ji ssoc -F'rem(CE t( ,v4 0( ,A'' s( )/ A a 



<ABCA1 



(5a) 



(5b) 



Comparison of the denominators in (4a) and (5a) show that in the 
former expression, an additional parameter Sydney enters; 
however, it is small compared to A:abcai (the mean values being 
2.42 and 95.18 respectively; see Table 5). In the numerator, the 
main quantitative difference between the two expressions is the 
remodeling flux, fc^j ssoc _Frern(CE 0 ,,^4 0 ,,Ay, versus the ApoA-I 
normalized flux, fc ( jj ssoc i ; rem(CE,,y4 OI ,A r oi)/^la. As shown in 
Figure 16, the latter has a flatter dependence on &ABCA1 as 
well as less variability due to other parameters. As a result, the 
ratio A/ p /A x allows for a more precise estimate of /:abcai 
compared to A/ p above. In conclusion, the analysis shows that the 



stronger inverse relationship shown in Figure 15 panel B can be 
attributed to the normalization of the remodeling flux by ApoA-I. 

The simulation results may further explain why, in some 
literature studies, % lipid-poor ApoA-I (note that the absolute 
concentration of lipid-poor ApoA-I can be experimentally 
estimated by assays that measure pre-/?! [61]) has been proposed 
as a risk factor, as well as how increased % lipid-poor ApoA-I 
could be associated with CVD risk. Our proposal of using the % 
lipid-poor ApoA-I as a surrogate measure for ABCA1 activity is in 
concordance with the previous suggestion by Asztalos et al [62] 
that the ratio pre-/?i/ai is a measure of the efficiency of RCT: a 
decrease in this ratio has been thought to reflect an enhanced 
RCT [62,63]. In addition, our finding of the inverse correlation 
between % lipid-poor ApoA-I and ABCA1 activity may explain 
the observation that increased fractional pre-/?i is associated with 
increased maximum intima-media thickness in both diabetics [64] 
and non-diabetic subjects [65], as well as being associated with an 
increased risk for coronary heart disease and myocardial 
infarctions [66]. 

Future directions 

We foresee a number of potential future applications of the 
LMK model in the context of drug discovery and development, 
including the following: 

1 . Confirmation of a molecule's mechanism of action: this can be 
done by checking the clinically observed changes in lipoprotein 
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Figure 4. The fit of the model to the calibration data: FCR of ApoA-l versus HDL-C/ApoA-l ratio. The data sources are: Brinton et al [43], 
Ikewaki et al [45], Schaefer et al [44]. The piecewise linear fit and the confidence interval are discussed in the Methods section. The model simulation 
values are indicated by asterisk symbols, for the nominal subject and the heterozygote, homozygote of CETP mutation. 
doi:10.1371/journal.pcbi.1003509.g004 



measures against the model predictions. This is an important 
task, since a molecule that increases HDL-C may do so by 
modulating the RCT pathway not only on its intended target but 
may also have off-target effects. As the model shows, the choice 
of mechanism in raising HDL-C could be crucially important for 
whether or not it brings about cardiovascular benefit. 

2. Determining the right dosage schedule for maximum cholesterol 
removal: the LMK model could help to integrate the pharma- 
cokinetics of a molecule with the dynamics of HDL metabolism. 

3. Evaluating combinations of target modulations: the LMK 
model could help to address the question of the potential 
synergism between targets in the RCT pathway. 



4. Development of personalized health care (PHC) strategy: 
simulations of the model to generate virtual populations could 
be used to address the question of which patient subpopulations 
are most likely to benefit from a given therapy and how those 
subjects might be selected using plasma-based diagonostic tests. 

The LMK model is focused on capturing the dynamics of ApoA-I 
and CE transfers. However, extensions of the model to incorporate 
ApoA-II dynamics as well as explicitly representing triglyceride and 
phospholipid metabolism would be important for describing the effects 
of other drug classes, including the PPAR-a and y agonists [67,68] or 
synthetic phospholipids [69,70]. These remain topics for further 
research. 
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Figure 5. Model validation: simulation of ABCA1 and ApoA-l mutations compared with literature data for HDL-C (panel A) and 
ApoA-l (panel B). For the simulation results, the mean and the 95% confidence intervals are plotted. The data sources are Asztalos et al [46] and 
Santos et al [47]; the mean ± SD are shown. The model simulations of the mutation cases were obtained by taking the parameter values for the 
nominal subject and set /tabcai and r? to 50% and 0% of the nominal values respectively. 
doi:1 0.1 371 /journal.pcbi.1 003509.g005 



Conclusions 

We have developed a novel, in-silico model of lipoprotein 
metabolism focused on the reverse cholesterol transport pathway. 



The model incorporates important concepts of HDL biology, 
including the regeneration of lipid-poor ApoA-I via a-HDL 
remodeling processes, and has been calibrated using literature data 
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Figure 6. Model validation: Simulation of tracer kinetic experiment with labelled ApoA-l compared to experimental data. The data 
are obtained by digitization of tracer kinetics measurements carried out in 4 subjects and shown in Figure 3 of Ikewaki et al [48]. The model 
simulation corresponds to the nominal subject. 
doi:10.1371/journal.pcbi.1003509.g006 



from a wide variety of sources. The model has been further 
validated by simulating scenarios not considered in the calibration 
process. These include its ability to reproduce the levels of HDL-C 
and ApoA-I in hetero- and homozygous subjects with either 
ABCA1 or ApoA-I mutation and the observed biphasic kinetics of 
ApoA-I seen in tracer kinetics studies. This provides an increased 
confidence in the LMK model predictions with respect to 
modulations of these important targets and in the model's ability 
to simulate time-dependent scenarios. 

In this paper, we have illustrated the applications of the LMK 
model in comparing the two target modulations, CETP inhibition 
and ABCA1 up-regulation. The results drawn from our model 
provide a possible explanation for the non-efficacy of dalcetrapib in 
the dal-OUTCOMES trial [7] as well as suggesting that ABCA1 is a 
target that would increase the RCT rate. The model provides 
predictions on the biomarker changes as a result of ABCA1 target 
modulation. Furthermore, computational experiments using a virtual 
population have shown why the % lipid-poor ApoA-I, rather than the 
absolute concentration of lipid-poor ApoA-I, is a better biomarker for 
assessing the in-vivo ABCA1 activity. By integrating mechanistic 
concepts and data, the model provides a way to quantitatively 
evaluate and explore hypotheses of lipoprotein metabolism. 

Methods 

Model derivation 

Mass balance considerations. The LMK model (Figure 1) 
explicitly represents the mass balance of ApoA-I and CE molecules 



in plasma, whereas the mass balance of FC and PL molecules is 
represented implicitly. The input of ApoA-I to plasma reflects its 
synthesis rate, while the elimination of ApoA-I results from the 
excretion of lipid-poor ApoA-I by the kidney and holo-uptake of a- 
HDL particles by the liver. The remodeling of HDL particles by 
particle fusion, CETP, SRB1 and other processes leads to the 
recycling of ApoA-I from a-HDL particles to lipid-poor ApoA-I. 
Recycling influences the kinetics of ApoA-I in plasma but does not 
affect its mass balance. The input of CE to plasma reflects the 
rapid esterification of FC molecules in the nascent discs as they are 
converted to nascent spheres plus a small amount of CE which 
enters plasma during VLDL synthesis. The rate at which CE 
molecules appear in the a-HDL pool (via the nascent sphere) is 
defined in the LMK model as the RCT input rate and is assumed 
to equal the rate at which FC molecules are loaded onto the 
nascent discs (equation (2)). Elimination of CE from plasma results 
from holo-uptake and SRB 1 -mediated uptake of CE from all 
lipoprotein species. The CETP-mediated transfer of CE between 
a-HDL, VLDL and LDL does not affect the overall mass balance 
in plasma. 

FC and PL molecules are present on the surfaces of all spherical 
lipoprotein particles in plasma as well as on the membranes of red 
blood cells (RBCs) and other cells that are in contact with the 
plasma. Based on Shen€s model of lipoprotein structure (Shen et al 
[71] and Mazer et al [25]), we assume that FC is in rapid 
equilibrium between all of these species and that the amount of FC 
present on each particle surface (at equilibrium) is dependent only 
on the surface curvature, that is, the radius of the hydrophobic 
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Figure 7. The distribution of RCT rate and HDL-C and their correlation in the simulated virtual population. By drawing the parameters 
of the model from an uncorrected, multivariate normal distribution, a set of 2000 virtual patients is generated and the model simulations of RCT rate 
and HDL-C are shown. The right-hand axis represents the hypothetical inverse relationship between RCT rate and CVD risk. 
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core of the particle (as represented in equation (10), below). The 
FC needed for the surface of the nascent spheres is assumed to be 
provided by the large pool of FC present in blood, including 
RBCs, and is largely replenished by the HDL remodeling 
processes. It can be shown that the rate at which FC is eliminated 
from plasma via holo-uptake of HDL particles is very small 
compared to the RCT rate (<4%) and is therefore negligible from 
the perspective of mass balance. Similar considerations apply to 
the mass balance of PL. 

Particle size. The size of spherical a-HDL particles is 
computed in the model as follows. From the pool size of CE in 
a-HDL and the particle concentration, the number of CE 



molecules per HDL particle is given by: 



particle _ CE„(f) 



'CE 



N a {t)m c 



(6) 



Since CE a is expressed as an equivalent mass of FC and nic is the 
molecular weight of FC, Wq^^ 6 is appropriately determined. 
With an assumed ratio of TG/ CE = 0.13 in the core of HDL 
particles [25], we sum the volumes occupied by CE and TG to 
obtain the total core volume and determine the core radius (rcore) 
from it. Finally, the surface thickness ? = 20.2 A is added to the 



PLOS Computational Biology | www.ploscompbiol.org 



13 



March 2014 | Volume 10 | Issue 3 | e1 003509 



A Model of Lipoprotein Metabolism and Kinetics 



A 



Correlation coefficient =-0.32, p-value < 0.01 

o 10 r 



5 8 



LU 

O 
i 

_i 
Q 
X 

O 

CD 
-t— < 

03 
i_ 

CD 
O 

c 

03 



03 0 

u 0 
O 



• * t * \ 

• • Mm • mm mm* m * .a 




20 40 60 80 100 120 
HDL-C (mg/dL) 




Figure 8. The distribution of the clearance of HDL-CE and HDL-C and their correlation in the simulated virtual population. By 

drawing the parameters of the model from an uncorrelated, multivariate normal distribution, a set of 2000 virtual patients is generated and the 

model simulations of HDL-CE clearance rate and HDL-C are shown. 

doi:10.1371/journal.pcbi.1003509.g008 



core radius, giving rise to following expression for the particle 
diameter, d (in A): 



3 particle 



rcore=[j-n^—x(V CE +0.UxV TG ) , (7) 



1/3 



d=2x(i core +20.2). 



(8) 



Note that the molecular volumes Vqj; and VfQ are defined in 
Table 3. 

Remodeling flux. In the derivation of the remodeling flux, 
we compute the excess (or deficit) of ApoA-I compared to that 
derived using the updated Shen's model [25]. Given the pool size 
for AJj) and the particle concentration N x (t), we compute the 
number of ApoA-I molecules per particle: 



particle _ A a (t) 
7 ApoA-I~ N x (t)m A 



(9) 



Using the expression for rcore given in (7), the number of ApoA-I 
molecules needed to cover the surface is derived in the following 
manner [25]: firstly, the number of free cholesterol is computed, 



nc- 



(10) 



((^core + 20.2) 3 -r£ ore ) x exp(-84.4/(r C ore +20.2) -6.09) 



Then, the number of phospholipid molecules needed to cover the 
remaining surface area of the core is computed: 



4nrZ, 



7 PL = 



n C A C 



(11) 



where the cross-sectional surface areas of cholesterol and 
phospholipid are ^4^ = 39.1 A 2 and ^4pj^ = 68.5 A 2 respectively. 
The number of amino acids needed to cover the hydrophobic area 
exposed at the outer surface layer of HDL particle is: 



PLOS Computational Biology | www.ploscompbiol.org 



14 



March 2014 | Volume 10 | Issue 3 | e1 003509 



A Model of Lipoprotein Metabolism and Kinetics 



200 
T5 150 

5 100 

i 

q 50 

X 

0 

c 

12 
B 11 
1 10 

e 9 

8 




0 50 100 

CETP level (% of nominal) 




0 50 100 

CETP level (% of nominal) 



B 



5 200 

CD 

E 



i 

< 
o 

CL 

< 



o 
E 

Q_ 
I 

_l 
Q 
X 



100 
0 

, 20 
15 
10 
5 



0 




0 50 100 

CETP level (% of nominal) 




0 50 100 

CETP level (% of nominal) 



10 



0 



0 50 100 

CETP level (% of nominal) 



?g 2 
i— 

*1>1 



0 



0 50 100 

CETP level (% of nominal) 



Figure 9. Model predictions for the dependence of HDL measures (HDL-C, panel A; ApoA-l, panel B; HDL size, panel C; HDL particle 
concentration, panel D; lipid-poor ApoA-l, panel E) and RCT (panel F) on the CETP level. The model simulation curves were obtained by 
decreasing the 3 parameters associated with CETP activity (/i'hv TP >^hl TP '^lh TP ) f rom 100% to 0% of those corresponding to the nominal subject. For 
each prediction, the mean and the 95% confidence intervals are plotted. 
doi:10.1371/journal.pcbi.1003509.g009 



1 AA = 



47i(r C ore +20.2) 2 -«pl^PL 



*AA 



where the cross-sectional area of an amino acid A 



AA : 



(12) 



■ 15.6 A 



Since there are 243 amino acids in ApoA-I, under the further 



assumption that the weight fraction of ApoA-I in the HDL 
proteome is 60% [25], the number of ApoA-I molecules needed to 
cover the surface of HDL is given by: 



„Shen 
^ApoA-I" 



W AA x 
243 ' 



(13) 
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Figure 10. Simulation of CETP inhibition on a virtual population with low HDL-C (£40 mg/dL). Each virtual patient selected for the 
treatment simulation had its rate constants associated with CETP activity (A'Sy^P ^CETP anc | ^CETP^ d ecreasec | to 20% of their original values. 
doi:10.1371/journal.pcbi.1003509.g010 



Finally, the discrepancy between the number of ApoA-I on the 
HDL (9) and the number needed from the Shen model (13), is the 
excess (or deficit) ApoA-I. Given the HDL particle concentration 
N a , the following is the concentration of ApoA-I on ot-HDL which 
is available to dissociate as the remodeling flux: 



calculation of size d is given by equation (8) and the division by 
10 accounts for the conversion from A to nm: 



*holoW = ^holo+ /f holo x ( To ~ 7nm) 



(15) 



i r rem(CE a (/)^ ct (0,A f a(0) = 

m A x TV, x I » ApoA _ j « ApoA _ T 



(14) 



HDL holo-particle uptake. Our model allows for the 
possibility of a linear size dependence of the HDL holo-particle 
uptake rate. However, no prior assumption is made regarding 
the size dependency; using the calibration data, the sign and 
magnitude of the linear dependence is determined. In particular, 
the rate of holo-uptake has the following form, where the 



Parameter priors 

In this section, prior estimates of model parameters are given, 
including references to the original literature and the rationale for 
the choice of prior and the level of uncertainty. In a manner 
similar to a previously proposed Bayesian approach [35], 
uncertainty is increased by a factor /ijriap = 2 in the following 



• quantities that are measured in-vitro and mapped to the in-vivo 
context; 
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Figure 11. Model predictions for the dependence of HDL measures (HDL-C, panel A; ApoA-l, panel B; HDL size, panel C; HDL 
particle concentration, panel D; lipid-poor ApoA-l, panel E) and RCT (panel F) on ABCA1 activity. The model simulation curves were 
obtained by increasing the parameter representing ABCA1 activity (A'abcai) from 100% to 300% of the nominal subject. For each prediction, the 
mean and the 95% confidence intervals are plotted. 
doi:10.1371/journal.pcbi.1003509.g011 



quantities that are measured in a population with mutation(s) 
and mapped to normal subjects; 

quantities that result from pooling data obtained using distinct 
experimental techniques/ assumptions. 



No explicit prior correlations are assumed. 

Synthesis rate of ApoA-I (rP). The kinetics of ApoA-I were 
measured in n = 20 (11 males, 9 females) healthy subjects, with 
mean HDL-C = 46 mg/dL and ApoA-I= 115 mg/dL [44]. The 
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Figure 1 2. Simulation of ABCA1 up-regulation on a virtual population with low HDL-C (<40 mg/dL). Each virtual patient seleted for the 

treatment simulation had its ABCA1 activity W'ABCAl' ' ncrease d by 100% of its initial value. 

doi:10.1371/journal.pcbi.1003509.g012 



ApoA-I synthesis rate has been estimated to be 27.44 + 5.29 mg/ 
dL/day (mean±SD). Hence, we take r l P =21 A A ± 1.18 mg/dL/ 
day (mean±SEM). 

Rate of kidney elimination (^kidney)* There have been a 
number of papers describing the measurement of the FCR of 
pre-/?! [22,72]. However, the quantification of pre-/?! can be a 
challenging task and a more direct assessment of the clearance 
rate of lipid-poor ApoA-I is estimated by the FCR of ApoA-I in 
Tangier patients. In the model representation of homozygous 
Tangier patients (^ABCAl = ")> the FCR of ApoA-I equals the 
kidney clearance of lipid-poor ApoA-I. The hypothesis that 
kidneys is responsible for a large fraction of ApoA-I clearance is 
supported by the study of Braschi et al done using rabbits [73], 
where it was estimated that the kidneys contribute around 70% 
of total ApoA-I clearance. The residence times (RT) of ApoA-I 
in Tangier disease patients have been measured to be 0.30,0.13 
day in [74] and 0.22 day in [75]. It is assumed that in Tangier 
patients, ^kjdney = ' /RT. Thus, Sidney has been estimated to 



be 5.19+ 1.30 pool/day (mean±SEM). Because of the assump- 
tion made in mapping ApoA-I clearance measured in the 
Tangier patients to the normal population, we apply the factor 
/i ma p = 2 to give Sidney = 5.19 + 2.60 pool/day (mean±"' 
SEM). 

Dissociation rate constant of labile ApoA-I (^dj ssoc ). ^ n 

[76], fluorescence resonance energy transfer spectroscopy was 
used to quantify the rate of ApoA-I exchange. In this in-vitro 
set-up using synthetic rHDL incubated with 5-molar excess of 
lipid-free ApoA-I, the exponential relaxation time (defined as 
the time by which 50% of the exchange has occurred) was 
inferred to be 0.94 hour. This gives rise to the estimate of 
^(jjgsoQ = 17.7 pool/day [76]. The ApoA-I found on a-HDL 
particles can be divided into a tightly-bound pool [77] and a 
labile pool [78]. The study of the dissociation of ApoA-I 
molecules from the labile pool is carried out in [78], where 
surface plasmon resonance was used to study the kinetics of 
ApoA-I interaction with HDL particles. A two-state binding 
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Figure 13. Comparison of CETP inhibition with ABCA1 up-regulation: changes in RCT rate (panel A) and biomarkers (ApoA-l, panel 
B; HDL size, panel C; HDL particle concentration, panel D; LDL-C, panel E) versus the rise in HDL-C. The nominal subject is taken as the 
baseline. The model simulation of CETP inhibition is compared with literature data of CETP inhibitors, Dalcetrapib [50], Torcetrapib [51] and 
Anacetrapib [52]. 

doi:10.1371/journal.pcbi.1003509.g013 
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Figure 14. Correlations between HDL-C and HDL-P (panel A), and HDL-C and HDL size (panel B) in a virtual population of 2000 
subjects. 

doi:10.1371/journal.pcbi.1003509.g014 



model was used to describe the association and dissociation 
reactions and the rate parameters were identified from 
the time-course data. It has been found that for the pool of 
ApoA-I molecules that are bound to HDL particles in a stable 
conformation, the half-time of dissociation is around 3 minutes, 
corresponding to ^(jj ssoc «330 pool/day. Computing the 
mean and SD of the two estimates of ^(jj ssoc and using 
rtiiap = 2 to account for the fact that these values were 
measured in-vitro, we obtain ^(jj ssoc = 1 74 + 3 1 2 pool/day 
(mean±SEM). 

Rate constant of the lipidation of lipid-poor ApoA-I via 
ABCA1 (^ABCAl) - ^he m °del assumes that the lipidation 
of ApoA-I is initiated by ABCA1, leading to the formation of 



nascent discs and subsequently to nascent spheres (via 
LCAT). While LCAT is crucial for the esterification of 
free cholesterol to cholesteryl ester, ABCA1 activity is 
assumed to be rate-limiting in the formation of a-HDL. In 
the model, the rate at which the concatenation of 
processes leading from lipid-poor ApoA-I to mature, a-HDL 
is described by the ABCA1 activity, ^ABCAL Based on the 
size exclusion chromatographic technique for separating HDL 
into subclasses, the rate constant in the conversion of lipid- 
poor ApoA-I to the a-HDL pool has been estimated to be 
96.24 + 42.99 pool/day (mean±SD, n = 6) [22]. Thus, the 
mean and SEM is given by A'abcai =96.24+ 17.55 pool/day 
(mean±SEM). 



Correlation coefficient =-0.23, p-value < 0.01 
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Figure 16. Comparison of two expressions involving remodeling flux: %issoc xFrem (panel A) and ^jissoc x^rem/A (panel B). The 

distributions of absolute and ApoA-l adjusted remodeling flux in the virtual population are plotted against k^Q^ . The simulations of the nominal 
subject with only the parameter ^ABCAl variec l are shown as solid lines. 
doi:1 0.1 371 /journal.pcbi.1 003509.g01 6 



Stoichiometry of FC to ApoA-I in nascent discs (y) . In the 

model, y denotes the stoichiometry (or molar ratio) of FC to ApoA- 
I in nascent discs. Due to the model assumption that FC in nascent 
discs are all esterified and result in the formation of nascent 
spheres, y also equals the stoichiometry of CE to ApoA-I entering 
the a-HDL pool. Given these two interpretations of the model 
parameter y, there are alternative ways to estimate it from 
literature data. In [79], plasma was fractionated using two- 
dimensional electrophoresis and the composition of pre-/?; was 
analyzed. The weight fraction of FC to ApoA-I was found to be 
0.16 + 0.04 (n = 4). Using the given molecular weights of ApoA- 
I and FC, we obtain the molar ratio of FC/ApoA-I = 1 1.7+ 1.5 
(mean±SEM). An alternative estimate of y is obtained using 
the estimates for the rate of cholesterol esterification to HDL- 
CE (77.2+12.5 mg/dL/day, n = 3) [20] as well as the 
production rate of a-HDL from pre-/?! (882 + 540 mg/dL/ 
day, n = 6) [22]. Thus, using these sets of data y is estimated to 
be 6.46 + 2.22 (mean ± SEM). Finally, an in-vitro experiment 
has been carried out to characterize the composition of nascent 
HDL (nHDL) formed by the action of ABC A 1 on ApoA-I [80]. 
For the small nHDL formed (diameter =7.5 nm), the particles 
were found to contain, on average, 2 ApoA-I and 9 total 
cholesterol. This gives the estimate for y of 4.5. Thus, 
combining all three estimates we get y = 7.55+ 1.97 (mean±"' 
SEM). Using the factor /%iap = 2 to take into account that in- 
vitro estimates were used, we obtain the final estimate of 
y = 7.55±3.94 (mean±SEM). 

Rate constant of CE transfer from HDL to VLDL 



(k c 



In [20], the unidirectional movement of CE from 



HDL to VLDL was quantified using a two -pool model for CE in 
the Apo B- 1 00 particle classes (VLDL and LDL) and a single pool 
for CE in HDL particles. Using the data from n = 3 subjects, 

^HV TP is estimated t0 be 2.24+1.31 pool/day (mean±SD). The 
rate of CE movement from HDL particles to VLDL has also been 
quantified in [21]. In this compartmental analysis, a 3-pool model 
has been used for CE in Apo-B particles (VLDL, IDL and LDL) 
and bidirectionality of transfer has been assumed between HDL 
and VLDL as well as between HDL and LDL. Due to the fact that 
our model does not account for IDL, the CE transfers for this 
density class are pooled into those of VLDL. For n = 7 normal 
subjects, the net transfers of CE from HDL to VLDL and IDL are 
36.5 + 9.5 mg/dL/day. Normalizing by the individual concentra- 
tions of HDL-CE, this gives the estimate of £^y TP = 1.14 + 0.37 
pool/day (mean±SD). By pooling the data sets from both 
Ouguerram et al [20] and Schwartz et al [21] and using 
^map = 2 to account for the difference in the structures of 

(""FTP 

compartmental models, we obtain: ^jjy =1.47 + 0.58 pool/ 
day (mean ± SEM). 

Rate constant of CE transfer from HDL to LDL 

(*hl TP )- In [20], the movement of CE from HDL to LDL was 
quantified using a two -pool model for CE in the Apo B-100 
particle classes (VLDL and LDL) and a single pool for CE in HDL 

CFTP 

particles. Using the data from the 3 subjects, * s esumate d 

to be 9.20 + 2.54 pool/day (mean±SD). A different estimate is 
obtained using the data for n = l normal subjects given in [21]: by 
dividing the rates of CE movement from HDL particles to LDL by 



the concentrations of HDL-CE at the individual level, ^j^l^ * s 
estimated to be 3.88 + 2.1 1 pool/day (mean±SD). By pooling the 
data sets from both Ouguerram et al [20] and Schwartz et al [21] 
and using ^map = 2 to account for the difference in the structures 

of compartmental models, we obtain: ^jjj^ =5.47 + 2.05 pool/ 
day (mean ± SEM). 

Rate constant of CE transfer from LDL to HDL 
(*LH TP )- In P°]> the rate of GE transfer from LDL to HDL 
was quantified using a two -pool model for CE in the Apo B-100 
particle classes (VLDL and LDL) and a single pool for CE in HDL 
particles. Using the data from the 3 subjects, ' s est i mate d 

to be 3.12+ 1.17 pool/day (mean±SD). An alternative estimate is 
obtained using the data for n = l normal subjects given in [21]: by 
dividing the rates of CE movement from LDL particles to HDL by 

the concentrations of LDL-CE at the individual level, ^LH * s 
estimated to be 1.49 + 0.64 pool/day (mean±SD). By pooling the 
data sets from both Ouguerram et al [20] and Schwartz et al [21] 
and using ^map = 2 to account for the difference in the structures 

of compartmental models, we obtain: = 1-98 + 0.70 pool/ 

day (mean ± SEM). 

Rate constant of transfer of CE from VLDL to LDL 

(ky L ). In [20], the rate constant of CE transfer from VLDL to 
LDL (due primarily to lipolysis) was inferred in « = 3 normal 
subjects: this gives rise to the parameter estimate &vl = 7.52 + 0.94 
pool/day (mean ± SEM). 

Flux of CE to VLDL (f-y LDL ) . In [20] , the flux of CE into the 
VLDL pool was inferred in « = 3 normal subjects to be 
0.96 + 0.8 mg/dL/day (mean±SD). The cholesteryl ester pro- 
duction to VLDL was also measured by Schwartz et al in [21], but 
due to the large uncertainty as represented by the greater than 
100% SD in some of the individual data, these values have not 

been used. Thus, we take ^ LDL = 0.96 + 0.46 mg/dL/day 
(mean±SEM). 

Rate constant of CE elimination from VLDL (&^ DL ). In 

[20], the rate constant of CE elimination from the VLDL pool was 
inferred in w = 3 normal subjects to be 0.88 + 0.64 pool/day 
(mean±SD). The quantification of CE elimination rate from 
VLDL was also carried out by Schwartz et al in [21], but the mean 
of the data was not shown in the paper because most values were 
undefined (fractional SD >80%). Thus, we use only the values 
given by Ouguerram et al [20] and take A^ u L t DL = 0.88 + 0.37 pool/ 
day (mean ± SEM). 

Rate constant of CE elimination from LDL (k^ L ). In 
[20], the rate constant of CE elimination from the LDL pool was 
inferred in w = 3 normal subjects to be 0.67 + 0.14 pool/day 
(mean±SD). The quantification of CE elimination rate from LDL 
was also carried out by Schwartz et al in [21], but the mean of the 
flux to the extra-hepatic pool was not shown in the paper because 
most values were undefined (fractional SD >80%). Thus, we use 
only the value given by Ouguerram et al [20] and take 
A;LDL =0 67±0 08 pool / day (mean±SEM). 

Rate constant of SRB1 -mediated CE elimination from 
HDL (£™Bi)- In [20], the rate constant of selective CE 
elimination from the HDL pool was inferred in w = 3 normal 
subjects to be 0.31 +0.20 pool/day (mean±SD). The selective 
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Table 6. Calibration data: HDL-C and ApoA-l in normal and CETP deficient subjects. 



Type 


Data 


Inazu [81] 
(mean±SD) 


Yamashita [82] 
(mean±SD) 


Asztalos [83] 
(mean±SD) 


Pooled (mean+SEM) 


Normal Subjects 




1 6 


20 


50 


86 




HDL-C 


52.9±13.9 


50±8 


52±14 


52±1 




ApoA-l 


124±21 


140.9±16.1 


144±29 


139±4 


Heterozygotes of CETP deficiency 


n 


20 


15 


5 


40 




HDL-C 


66±15 


84±25 


85 ±26 


75±4 




ApoA-l 


149±43 


155.3±22.1 


154±25 


152±5 


Homozygotes of CETP deficiency 


n 


10 


4 


9 


23 




HDL-C 


163.7±39 


193±28 


157±29 


166±7 




ApoA-l 


213±47 


233.5±22.3 


252±25 


232±8 



All concentrations are given in mg/dL. 
doi:1 0.1 371 /journal.pcbi.1 003509.t006 



uptake of HDL-CE by the liver was not observed in Schwartz 
et al in [21]. Hence, we use only the value given by 

Ouguerram et al [20] and take ^§^g| =0.31+0.12 pool/day 
(mean + SEM). 

Rate constant of size-independent holo-particle uptake 
for h-HDL (Af j ). The FCR of ApoA-I in the 0£-HDL pool 
has been estimated to be 0.112 + 0.026 pool/day (mean±SD, 
n = 6) using HDL subclasses separated with size exclusion 
chromatography [22]. In another reference [72], using a 
separation technique based on agarose gel electrophoresis, the 
FCR of ApoA-I in the a-HDL pool has been estimated to be 
0.17 + 0.039 pool/day (mean±SD, n — 6). Thus, by pooling the 
data we obtain fc£ j =0.14 + 0.013 pool/day (mean±SEM). 

Finally, using the factor /%mp = 2, we get ^olo = 0.14 + 0.026 
pool/day (mean±SEM). 

Calibration data 

In this section, we give the quantitative values and references for 
the data used in the calibration procedure. 

CETP mutation. Lipoprotein data for CETP mutation 
patients were taken from 3 sources: Inazu et al [81], Yamashita 
et al [82] and Asztalos et al [83] and were pooled to yield the mean 
and SEM. In particular, both HDL-C and ApoA-I were available 
for each of the 3 data sources shown in Table 6. The values for 
LDL-CE and VLDL-CE for the control subjects were not given in 
the references for CETP mutation [81-83]; hence, the values from 
Ouguerram et al [20] and Schwartz et al [2 1] were used in Table 7. 
The value of LDL-CE for CETP heterozygotes was also used for 
calibrating the model, which was estimated using the given value 



Table 7. Calibration data: CE in ApoB particles. 





Type 


Data [20,21] 


Pooled (mean+SEM) 


Normal subjects 


n 


10 




LDL-CE 


83±3 mg/dL 




VLDL-CE 


5±1 mg/dL 


Heterozygotes of CETP 
deficiency 


n 


23 




LDL-CE (LDL-C x0.7) 


72 ±9 mg/dL 



of LDL-C with an assumption on the ratio of FC/CE for LDL 
particles. In Asztalos et al [83], the ratio of FC/CE for all 
lipoprotein particle classes of CETP heterozygotes was given as 
0.37, which gives CE = 0.73xTC. Given the approximations 
made, we assumed that LDL-CE = LDL-C x (0.7 + 0.05). There 
is literature data for LDL-CE and/or VLDL-CE in CETP 
homozygotes [39-41], but these were not used in the calibration 
process due to the known inconsistency of the model in lacking the 
/J-LCAT activity [42]. 

Cholesteryl ester flux. There are 2 literature data sources 
on the in-vivo flux of CE between HDL and ApoB-containing 
particles: Ouguerram et al [20] and Schwartz et al [21]. While both 
data sets have been used in estimating the prior distribution of 
parameters involved in the exchanges of CE (see Methods section), 
for the purpose of model calibration a choice between the 2 
disparate data needed to be made. Given that the model structure 
in the description of CE exchanges between HDL, LDL and 
VLDL is based on that of Ouguerram et al [20] , it was decided that 
the same calibration data should be taken for the CE fluxes. The 
values used are shown in Table 8. 

Fractional catabolic rate of apoA-I. Brinton et al have 
shown that a strong association exists between the fractional 
catabolic rate (FCR) of ApoA-I and the estimated HDL size, using 
a surrogate marker [43] (Brinton et al used HDL-C /(ApoA- 
I+ApoA-II), we re-analyzed their data taking HDL-C /ApoA-I as 
the surrogate marker). Their work has demonstrated that as much 
as 70% of the variability in the FCR of ApoA-I may be attributed 

HDL-C 

to variations in HDL size, as estimated using [431. This 

°ApoA-I L J 

finding is corroborated with the individual data of ApoA-I 
metabolism from Schaefer et al [44] measured in healthy 
volunteers. Finally, FCR of ApoA-I was also studied by Ikewaki 
et al [45] in comparing CETP mutation subjects with controls. All 



Table 8. Calibration data: CE fluxes. 



doi:10.1371/journal.pcbi.1003509.t007 







CE flux (mean+SEM) 


Data (mg/dL/day) 


HDL to VLDL 


64±2 


HDL to LDL 


268±51 


LDL to HDL 


266 ±52 



doi:1 0.1 371 /journal.pcbi.1 003509.t008 
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these data consistently show that the ApoA-I FCR exhibits an 
inverse relationship with the surrogate measure of HDL size. 
However, the trend does not continue indefinitely: even for CETP 
homozygous subjects with very large particles, ApoA-I FCR 
appears to reach a minimum of 0.135 pool/day. Hence, we have 
tried to describe the data in the simplest way, assuming a linear 

HDL — C 

relationship between ApoA-I FCR and — — — -, with a lower 

bound of 0.135. Using this linear assumption, the fit and 
confidence interval was computed and shown in Figure 4. The 
piecewise linear fit to the combined data is given by 

/HDL-CA 

FCR ApoA _T=max(0.50-0.63 x ( apoA-l )' 0 - 135 ^ the 
mean confidence interval corresponding to 1 SD in the line fit is 
approximately ±0.065 pool/day. 

Model calibration 

In this work, we assume that both the parameter prior and the 
data error are normally distributed. We employ the methodology 
of maximum a posteriori (MAP) [31] to combine the prior information 
with calibration data. Due to the conjugacy property [31] of the 
distributions, the posterior also has a normal distribution and the 
MAP solution is obtained by solving a nonlinear least squares 
problem. In our model, most of the parameters have an 
informative prior. For the set of parameters for which an 
informative prior is available, let k pr j or denote the expected 
value of the prior distribution; otherwise, set k pr j or = 0 to 

represent the lack of information. We take the covariance matrix 
for the prior distribution CY. to have a diagonal structure: for 
parameters kj that have an informative prior, (C^)^- is the 
variance of the prior distribution; for parameters that have an 
uninformative prior, (C^.)^- = go. That is, the prior distribution is 
assumed to be of the form [33]: 



p(k)cj:exp(^-^(k-k prior ) T C k \k-k prior )y (16) 

Let deR. m denote the vector of calibration data and G(k) the 
nonlinear mapping from model parameters to the observation, 
representing the model simulation of the data. Let denote the 
covariance matrix for the data. Hence, the conditional distribution 
of the data given the model parameter k is [33] : 

f(d\k) ocexp ( - l - (G(k) -d) T C d \G(k) -d^j. (17) 

Thus, the posterior distribution q{k\d) for the model parameters is 
given by 



the MAP solution is the minimizer: 



q(k\d) rjzf(d\k)p(k) = exp ^ - - (G(k) - d) T Cj 1 (G(k) - d) - 

1 

2 



^ (k kp r i or ) Cf, {k kprior) 



(18) 



To find the MAP solution, the following nonlinear least squares 
problem is solved: with the objective function defined as, 



X 2 (k) = 

(G(k) -d) T C d l (G(k) -d) + (k- k prior ) T C k 1 (k - k prior ), 



(19) 



C MAP* 



minx (k)- 
k 



(20) 



Using parameter priors as given in Table 5 and calibration data as 
described in the previous section, the nonlinear least-squares 
problem was solved using genetic algorithm ga from the Matlab® 
Global Optimization Toolbox of MathWorks (http://www. 
mathworks.com/) to obtain ^map- In particular, the hybrid option 
was selected: 1 00 generations of the genetic algorithm was run 
with a PopulationSize = 500, followed by constrained minimiza- 
tion (fmincon) using the setting MaxFunEvals = 1 0000, Maxl- 
ter= 1000. In all numerical integration of ODEs, the relative and 
absolute tolerances were set to 10 9 . 

Estimation of 95% confidence intervals 

The confidence interval is estimated using the following 
procedure: parameters are sampled around ^map and for each 
parameter the A% 2 (with respect to its minimum, X 2 (^MAp)) * s 
computed according to the expression (19). An estimate of the 
confidence region is obtained by examining the set of all parameters 
that lie within A/ 2 < <5, where S is computed from the number of 
degrees of freedom (df) and the desired confidence level [84]. Using 
df=29 for the model and choosing the 95% confidence level, 
<5 = 42.557. A set of 1000 parameters satisfying A# 2 <42.557 are 
selected in estimating confidence intervals shown in the paper. 

Simulation of tracer kinetic studies 

The model simulations of the tracer kinetic experiment with 
labelled ApoA-I and the calculation of the FCR of ApoA-I were 
carried out using the technique of complex variable differentiation 
[85]. In particular, a small quantity of imaginary number 
representing the radio-labelled dose of ApoA-I is added to the 
lip id-poor pool at the start of tracer experiment and the imaginary 
component of the numerical solution is extracted to represent 
the dose remaining in the two pools of ApoA-I (lipid-poor and 
a-HDL). This method relies on the complex extension of analytic 
functions from the real line, which can be easily implemented on 
the Matlab platform [85]. As compared to the finite-differencing 
approach, the complex variable methodology does not suffer from 
subtractive cancellation error and hence is more accurate [85]. 
While this approach has not been applied to tracer kinetic 
simulations, it has been applied to the sensitivity analysis of 
biological models [86,87]. 

Supporting Information 

File SI LMK model implementation in SimBiology and 
Matlab formats. 
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